n=15
alpha=1
beta=1
library(reshape2)
library(plotly)
## Loading required package: ggplot2
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(ggplot2)

#Cut into bins the data, but name each bin as the mean(left,right) of each bin.
cut_jp = function(uncut_vector,n_bins){
  labs = cut(uncut_vector,n_bins)
  cut_vec = data.frame(lower = as.numeric( sub("\\((.+),.*", "\\1", labs) ),
                         upper = as.numeric( sub("[^,]*,([^]]*)\\]", "\\1", labs) ))
  cut_vec$mean= rowMeans(cut_vec)
  cut_vec$bin_ix = cut(uncut_vector,n_bins, labels = FALSE)
  
  labs = levels(cut(uncut_vector, n_bins))
  bins_info = data.frame(lower = as.numeric( sub("\\((.+),.*", "\\1", labs) ),
                       upper = as.numeric( sub("[^,]*,([^]]*)\\]", "\\1", labs) ))
  bins_info$mean= rowMeans(bins_info)
  bins_info$ix=as.numeric(rownames(bins_info))
  
  output= list(cut_vector= cut_vec, bins_info=bins_info)
  #cut_vec= bin_names$new
  #cut_vec= bin_names$lower
  return(output)
}

1 Expected Number of Solutions

Expected number of solutions for a knapsack problem with 15 items. And weights and values are sampled from Dirichlet distributions:

\[ (v_i,...,v_n) \sim Dir(\alpha, ... , \alpha ) \] \[ (w_i,...,w_n) \sim Dir(\beta, ... , \beta ) \] with \(\alpha = 1\) and \(\beta=1\).

#Expected Number of Solutions (Analytical Calculation) 
grid_size=100

Np=seq(0,grid_size)
Nc=seq(0,grid_size)

#Calculates the expeted values and stores them in the prob array
prob=array(-1,dim=c(length(Nc),length(Np)))
for(npi in Np){
  for(nci in Nc){
    suma=0
    for(s in seq(1,n)){
      suma = suma + choose(n,s)*(1-pbeta(npi/grid_size, alpha*s, alpha*(n-s)))*pbeta(nci/grid_size,beta*s,beta*(n-s))
    }
    prob[nci+1,npi+1]=suma
  }
}
#Reshapes the Array and assigns an nprofit and ncapacity to each expected value.
eSol_data=melt(prob)
names(eSol_data)=c('nCapacity','nProfit','ESol')
eSol_data$nCapacity=eSol_data$nCapacity/grid_size
eSol_data$nProfit=eSol_data$nProfit/grid_size
#Plots the expected number of solutions for given paramaters: n, alpha and beta
p <- plot_ly(eSol_data,x=~nCapacity,y=~nProfit,z=~ESol, marker = list(color=~ESol, size=1) )%>%
  add_markers() %>%
  layout(scene = list(xaxis = list(title = 'nCapacity'),
                      yaxis = list(title = 'nProfit'),
                      zaxis = list(title = 'E(# of Solutions)')))
p

2 Constrainedness of Search (Gent, 96): Kappa

##Calculate Kappa based on the analytical results
bins = seq(0,100)/100

#eSol_data$kappa = 1 - log(eSol_data$ESol,base =2)/log(2^n,base=2)
eSol_data$kappa = 1 - log(eSol_data$ESol,base =2)/n

eSol_data$nCapacity_bin = cut(eSol_data$nCapacity,breaks =bins,labels=FALSE)
eSol_data$nProfit_bin = cut(eSol_data$nProfit,breaks =bins,labels=FALSE)

eSol_data$lnCP= log( (eSol_data$nCapacity_bin/100) / (eSol_data$nProfit_bin/100) )

\(\kappa\) kappa plotted in Normalised capacity and profit space

p <- plot_ly(eSol_data,x=~nCapacity,y=~nProfit,z=~kappa, marker = list(color=~kappa, size=1) )%>%
  add_markers() %>%
  layout(scene = list(xaxis = list(title = 'nCapacity'),
                      yaxis = list(title = 'nProfit'),
                      zaxis = list(title = 'Kappa')))
p
dataInput = eSol_data

plo = ggplot(dataInput, aes(x=nCapacity, y=nProfit,z=kappa))+
  geom_contour(aes(colour=stat(level)), bins=20)+
  scale_colour_gradient(low = "lightgoldenrod2", high = "red")+
  #guides(color = guide_legend(title = "Solvability probability"))+
  theme(plot.title = element_text(hjust = 0.5))+
  #ggtitle(paste0("Kappa parameter (n=",n,")"))+
  theme_light() +
  labs(x="Normalised Capacity",y="Normalised Profit")+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        plot.title = element_text(hjust = 0.5),strip.text.x = element_blank(),
        text = element_text(size=20))

plo
## Warning: Removed 201 rows containing non-finite values (stat_contour).

2.1 Kappa and Solvability (Phase Transition)

Solvability based on actual simulations and solver results based on Gamma(1,1) sampling over the weights and values and the previous analytical results for n=15.

#Get the standard Phase transition (Solvability) from the simulations
pt_standard = read.csv("Simulations Output/gamma_sim/simulation_data_n15_gamma1_1.csv", stringsAsFactors = FALSE)
p <- plot_ly(pt_standard,x=~nCapacity,y=~nProfit,z=~solvability, marker = list(color=~solvability, size=1) )%>%
  add_markers() %>%
  layout(scene = list(xaxis = list(title = 'nCapacity'),
                      yaxis = list(title = 'nProfit'),
                      zaxis = list(title = 'Solvability')))
p
#Join Simulations and kappa 
join_sim_ana = merge(pt_standard, eSol_data, by=c("nCapacity","nProfit"))

cut_result=cut_jp(join_sim_ana$kappa,100)
join_sim_ana$kappa_bin = cut_result$cut_vector$mean

join_sim_ana_bined= join_sim_ana %>% group_by(kappa_bin) %>% summarise(solvability = mean(solvability))
plo = ggplot(join_sim_ana_bined, aes(x=kappa_bin, y=solvability))+
  geom_point() +
  theme_light() +
  labs(x="Kappa",y="Solvability")+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        plot.title = element_text(hjust = 0.5),strip.text.x = element_blank(),
        text = element_text(size=20))
  

plo